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We reveal that local interactions in graphene allow novel spin liquids between the semi-metal 
and antiferromagnetic Mott insulating phases, identified with algebraic spin liquid and Z2 spin 
liquid, respectively. We argue that the algebraic spin liquid can be regarded as the two dimensional 
realization of one dimensional spin dynamics, where antiferromagnetic correlations show exactly the 
same power-law dependence as valence bond correlations. Nature of the Z2 spin liquid turns out 
to be d + id' singlet pairing, but time reversal symmetry is preserved, taking d + id' in one valley 
and d — id' in the other valley. We propose the quantized thermal valley Hall effect as an essential 
feature of this gapped spin liquid state. Quantum phase transitions among the semi-metal, algebraic 
spin liquid, and Z2 spin liquid are shown to be continuous while the transition from the Z2 spin 
liquid to the antiferromagnetic Mott insulator turns out to be the first order. We emphasize that 
both algebraic spin liquid and d ± id' Z2 spin liquid can be verified by the quantum Monte Carlo 
simulation, showing the enhanced symmetry in the algebraic spin liquid and the quantized thermal 
valley Hall effect in the Z2 spin liquid. 
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I. INTRODUCTION 

Interplay between the topological band structure and 
interaction drives one direction of research in modern 
condensed matter physics?^ where emergence of Dirac 
fcrmions is at the heart of the interplay. The original ex- 
ample is the system of one dimensional interacting elec- 
trons, where interactions become enhanced at low en- 
ergies, combined with one dimensionality, and electron 
fractionalization results, giving rise to a new state of 
matter, dubbed as the Tomonaga-Luttinger liquid^ An 
interesting aspect is that such fractionalized excitations 
as spinons and holons are identified with topological ex- 
citations, carrying fermion quantum numbers associated 
with the topological structure of the Dirac theory^ 

A recent study based on the quantum Monte Carlo 
simulation 5 has argued that essentially the same phe- 
nomenon as electron fractionalization in the Tomonaga- 
Luttinger liquid may happen in two dimensions when lo- 
cal interactions arc introduced in the graphene structure. 
This study claimed existence of a paramagnetic Mott in- 
sulator with a spin gap between the semi-metal and anti- 
ferromagnetic Mott insulating phases. Immediately, the 
nature of the spin gapped Mott state has been suggested 
to be an s-wave spin-singlet pairing order between next 
nearest neighbor spins* 6 - - — thus identified with a Z2 spin 
liquid. We point out other scenario a 9 ' 10 for the nature of 
the spin liquid state. 

In the present study we revisit this problem, the na- 
ture of possible spin liquids in the Hubbard model on 
the graphene structure. An important point of our study 
is to solve the Hubbard model directly beyond recent 
studies* 6 - - — where an additional energy scale was intro- 
duced to describe the spin-singlet pairing order. The 
SU(2) slave-rotor representation, invented by one of the 
authorSfii is at the heart of the methodology, where ex- 



change correlations via virtual processes are naturally 
caught to allow spin singlet-pairing. One may regard 
the SU(2) slave-rotor theory of the Hubbard model as an 
analogue of the SU(2) slave-boson theory^ 2 - for the t-J 
model. 

We find two kinds of spin liquids, identified with an 
algebraic spin liquid and a Z2 spin liquid, respectively, 
between the semi- metal and antiferromagnetic phases. 
We argue that the algebraic spin liquicU^^ can be re- 
garded as the two dimensional realization of one dimen- 
sional spin dynamics, where antiferromagnetic correla- 
tions show exactly the same power-law dependence as va- 
lence bond correlations ! 15 ' 16 Increasing interactions, pair- 
ing correlations between nearest neighbor sites become 
enhanced. As a result, the algebraic spin liquid is shown 
to turn into a gapped spin liquid state, where the spin gap 
results from d + id' singlet pairing, believed to originate 
from the interplay between the topological structure and 
interaction. Actually, this pairing symmetry solution has 
been argued to be stable, based on an effective model in 
the weak coupling approach. 17,18 However, we argue that 
time reversal symmetry is preserved, taking the d + id' 
singlet pairing to one valley while the d — id' pairing to 
another. This assignment turns out to be essential in or- 
der to have a fully gapped spectrum because the d + id' 
singlet pairing order parameter in one valley vanishes in 
the other valley, allowing the gapless Dirac spectrum. We 
propose the quantized thermal valley Hall effec t 19 ' 20 for 
the fingerprint of this gapped Z 2 spin liquid. We dis- 
cuss the nature of quantum phase transitions beyond the 
mean-field analysis. 

We would like to emphasize that appearance of both 
algebraic spin liquid and d ± id' Z2 spin liquid can be 
verified by the quantum Monte Carlo simulation in prin- 
ciple. The fingerprint of the algebraic spin liquid is an 
enhanced symmetry, giving rise to the same power-law 
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dependence between antiferromagnetic and valence bond 
correlations. The hallmark of the d ± id' Z 2 spin liquid 
is the quantized thermal valley Hall effect, as mentioned 
above. We hope that the present study motivates quan- 
tum Monte Carlo simulation researchers to calculate such 
correlation functions. 

The present paper is organized as follows. In Sec. II 
we present the SU(2) slave-rotor theory of the Hubbard 
model. General mean field equations are also obtained. 
The mean field analysis of possible phase transitions is 
presented in Sec. III. In Sec. IV summary and discussion 
are presented. 



II. SU(2) SLAVE-ROTOR THEORY OF THE 
HUBBARD MODEL 



A. Formulation 

We start from the Hubbard model on the honeycomb 
lattice 

H = -t ^ 4a C ja + H ' C - + U E ni T n 4> (!) 

(ij}<7 I 

where c; CT {c\ a ) is the annihilation (creation) operator 
for an electron at site i with spin a. t is the hopping 
integral, and U is the on-site Coulomb interaction, where 
Uia = cla-Cia represents the density of electrons with spin 
a. 

Introducing the Nambu-spinor representation 
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and performing the Hubbard-Stratonovich transforma- 
tion for the pairing, density (singlet) and magnetic 
(triplet) interaction channels, we obtain an effective La- 
grangian 

C = y"^t(drl - ^(J z )ipz - ty^^l(T z ipj +H.c. 
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Here, $ t and <fi are associated with pairing- 
fluctuation and density-excitation potentials, introduced 
to decouple the charge channel, m, is an effective mag- 
netic field, which decouples the spin channel. n c and k s 
are introduced for decoupling between singlet and triplet 
interactions in respect that we do not know how they be- 
come renormalized at low energies. One may regard these 



two decoupling coefficients as phenomenological param- 
eters to overcome the mean-field approximation. Several 
examples for decoupling are shown in appendix A. We 
emphasize that both the semi-metal to algebraic spin liq- 
uid and the algebraic spin liquid to Z 2 spin liquid phase 
transitions are shown not to depend on such phenomeno- 
logical parameters, where both the critical points are 
determined with the combination between U/t and k c . 
Only the Z 2 spin liquid to antiferromagnetic transition 
turns out to depend on such parameters. We should be 
careful in determining this phase transition, comparing 
various cases (appendix A) with each other. 

The SU(2) slave-rotor representation 11 means to write 
down an electron field as a composite field in terms of a 
charge-neutral spinon field and a spinless holon field 
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where Fi = ( 



a fermion operator in the Nambu 
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representation, and Zi is an SU(2) matrix 
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Here, zi a is a boson operator, satisfying the unimodular 
(rotor) constraint, zLz^ + 2 J, Zn = 1. 

The key point of the slave-rotor representation^! is 
to extract out collective charge dynamics explicitly from 
correlated electrons. Such charge fluctuations are identi- 
fied with zero sound modes in the case of short range in- 
teractions while plasmon modes in the case of long range 
interactions. Actually, one can check that the dispersion 
of the rotor variable (Zif) is exactly the same as that of 
such collective charge excitations. 

In the slave-rotor theory the Mott transition is de- 
scribed by gapping of rotor excitations^ Until now, the 
Mott transition has not been achieved successfully, based 
on the diagrammatic approach starting from the Fermi 
liquid theory. In this respect the slave-rotor theory can 
be regarded as an effective field theory, not bad for the 
Mott transition. 

Resorting to the SU(2) slave rotor representation in 
Eq. ©, we rewrite the effective Lagrangian Eq. (2) as 
follows 
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FIG. 1: Graphene. A distance between two next nearest 
neighbor sites is chosen as the length unit, ai and a.2 are 
primitive translation vectors. 81, 82, and #3 are three nearest 
neighbor bonds. 

It is not difficult to see equivalence between the SU(2) 
slave-rotor effective Lagrangian and the decoupled Hub- 
bard model [Eq. (2)]. Integrating over field variables of 
Xij and Yij, and shifting fi^-cr as 

fti-a+iZid T Zj - ifiZi<T z Zj, 

where fli = <&f , ifii) is the pseudospin potential field, 
we recover Eq. (2) exactly with an introduction of an 
electron field ZjFi — > ipi. This procedure is well de- 
scribed in the previous study.™ An important feature in 
the SU(2) slave-rotor description is appearance of pair- 
ing correlations between nearest neighbor electrons, given 
by off diagonal hopping in which results from on-site 
pairing (virtual) fluctuations, captured by the off diago- 



nal variable zn of the SU(2) matrix field Z;. We note that 
the diagonal rotor field corresponds to the zero sound 
mode, giving rise to the Mott transition via gapping of 
their fluctuations. The additional boson rotor variable 
Zii allows us to catch super-exchange correlations in the 
Mott transition. But, the appearance of pairing correla- 
tions does not necessarily lead to superconductivity be- 
cause their global coherence, described by condensation 
of SU(2) matrix holons, is not guaranteed. The similar 
situation happens in the SU(2) slave-boson theory^ of 
the t-J model. 

B. Mean- field ansatz 

We perform the mean-field analysis, taking the follow- 
ing ansatz 

where 8 denotes the bond between the nearest neighbor 
sites. In the honeycomb lattice there are three near- 
est neighbor bonds. We choose wg = wjs-i vg — vQg, 
wg = wjg, and vg = v($, where 75 and (g are symmetric 
factors for the hopping parameter w (w) and the pairing 
order parameter v (v). The choice for 75 and Cs depends 
on the symmetry of the considered phase. For example, 
the s-wave pairing symmetry is given by Q = (1, 1, 1), the 
d x 2_ y 2-w&ve symmetry £5 = (— i,— and the d xy - 
wave symmetry = {— h,— h,0). For the magnetic or- 
der parameter rrii we choose an antiferromagnetic ansatz 
rrii = (— l) l m. 

Particle-hole symmetry at half filling results in /i + 
i(fi = while pairing potentials of $f and $| vanish in 
the mean-field level. Then, we obtain a general expres- 
sion for the free energy 
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where 7(k) = ^ 5 75 exp(ir,5 ■ k) is the energy dispersion for spinons and holons, and £(k) = J2s C<5 exp(ira • k) is 
associated with the pairing potential. Y)g is performed in the unit cell. A is a Lagrange multiplier field, introduced 
to keep the slave-rotor constraint. N is the total number of sites. 

Minimizing the free energy, we obtain fully self-consistent equations for order parameters 
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In this study our objective is to reveal the phase structure of the Hubbard model on the honeycomb lattice. It 
is convenient to take the zero temperature limit. Performing the Matsubara frequency summation, we obtain self- 
consistent mean-held equations at zero temperature 



w ^ | 7 (k)| 2 , w ^ |T(k)l = 
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where 

E{k) = ^ 2 |7(k)| 2 +« 2 |C(k)| 2 , (24) 
D(k, to) = y w 2 | 7 (k)| 2 + («|C(k)| + ^) 2 (25) 

are holon and spinon energy spectra in the presence of 
pairing and antiferromagnetism, respectively. 

Considering symmetry, it is natural to take into ac- 



count spatially uniform hopping 

1 /3 

| 7 (k)| 2 = 3 + 2cos(£g+4cos(-fc y )cos(^-fc ;c ). (26) 

On the other hand, the s-wave pairing potential is not 
allowed due to repulsive interactions. Counting the lat- 
tice symmetry of the honeycomb structure, the next can- 
didate will be d x 2_ y 2 or d xy for nearest neighbor singlet 
pairing 1 ^. We introduce a general combination of d x i_ y i- 
and d^y-wave pairing for the pairing term C(k) 

|C(k)| 2 = I cos^C^ (k) -Msin^C^k)! 2 , (27) 
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where 9 is a combination factor, and Ca; 2 -y 2 [Cxy) ls the 
d x 2_ y 2 (d xy ) -wave symmetry function 

(x2-y2(k x ,k v ) = e~ lh n - cos(^) 

(xy{k x ,k v ) = ie^sin(y). 

For = ±7r/3 this pairing symmetry becomes d ± id'. 
We also consider the d + d'-wave pairing symmetry 

|C(k)| 2 = Icos^C^k) +sm(9)\<; xy (k)\ 2 , (28) 

but this pairing order turns out to be not a solution of the 
mean-field equations. If one tunes k c and k s parameters, 
he can make this pairing symmetry a solution. However, 
this solution does not give the lowest free energy, com- 
pared with the d + id' pairing solution, consistent with 
earlier studiesj 17 i 18 

One may criticize the ansatz for uniform hopping in 
this paper because such an assumption excludes possible 
dimerized phases a priori. Actually, the J\ — J2 Heisen- 
berg model 

H = J l &i ' &3 + ^2 ^ Sk ■ Si 

(ij) {(«» 

has shown several types of dimerized phases when the 
ratio of J2/J1 is beyond a certain critical value^^, ap- 
proximately given by J2/J1 ~ 0.2 ~ 0.3. Here, the first 
term represents the exchange interaction between near- 
est neighbor spins, and the second expresses that between 
next nearest neighbor ones. This model can be derived 
from the Hubbard model, resorting to the degenerate per- 
turbation theory in the t/U — > limit 23 , where each pa- 
rameter is given by£ 

^«{i-4(i) S }, A = 4( (i) 3 

up to the fourth order process. Then, the J2/J1 ratio 
can be expressed in terms of U/t as follows 

•h = 1 

j a (u/tr-4- 

It was argued that higher order terms such as third 
neighbor and ring exchange terms may be ignored be- 
cause third neighbor exchange terms are not frustrat- 
ing, just renormalizing the J\ term effectively, while the 
ring exchange term is expected to be small^ However, 
the role of the ring exchange term has been also studied 
carefully . 24 ' 25 

An antiferromagnetic phase has been reported in 
J2/J1 < (J 2 /Ji)af « 0.08£2£ This corresponds to 
{U/t)AF ~ 4.3, consistent with the result of the quan- 
tum Monte Carlo simulation. 5 Increasing frustration, the 



antiferromagnetic order disappears, and a paramagnetic 
Mott insulating state results, identified with a certain 
type of Z2 spin liquids. Such a spin-gapped state turns 
out to evolve into a dimerized phase with either trans- 
lational or rotational symmetry breaking near J2/ J\ ~ 
0.2 ~ 0.3j&22 It was reported that the spin liquid state 
turns into a dimerized phase with three-fold degeneracy 
around J2/J1 ~ 0.3, where it breaks the C3 symme- 
try but preserving the translational symmetry^ On the 
other hand, the plaquette order was claimed to appear 
near J2/J1 ~ 0.2 before the dimerized phase, breaking 
the translational symmetry onlyj22 An important point 
is that if we translate the critical J2/J1 value in terms of 
U/t of the Hubbard model, J2/J1 ~ 0.3 corresponds to 
U/t fa 2.7 and J2/J1 « 0.2, U/t ps 3.0. Comparing these 
critical values with the critical value for the semi-metal 
to spin liquid transition in the quantum Monte Carlo 
simulation-, the Mott critical value given by U/t ps 3.5 
turns out to be larger than those for dimerized phases. 
This means that the semi-metal phase will appear before 
reaching such dimerized phases in the Hubbard model 
owing to charge fluctuations, not introduced in the J\ — J 2 
Heisenberg model. In other words, the J\ — J 2 model 
seems to be an effective low energy model only in the 
limit of U/t — > 00 while physics of such a model will be 
different from that of the Hubbard model in the small 
U ft case. 

However, it should be pointed out that these critical 
values cannot be guaranteed. Thus, we cannot exclude 
the possibility of dimerization near the Mott criticality of 
the Hubbard model completely. In addition, introduction 
of the next nearest neighbor hopping t' will favor dimer- 
ization. In this respect it will be the best interpretation 
that the spin liquid physics may appear at finite tem- 
peratures at least, actually observed from the quantum 
Monte Carlo simulation. 



III. SADDLE-POINT ANALYSIS 
A. From semi-metal to algebraic spin liquid 

The semi-metal phase is described by condensation of 
holons (z a ) 7^ with v = v = m = 0. Considering the 
symmetry factor 75 = (1, 1, 1), the condensation occurs 
when the effective chemical potential given by the La- 
grange multiplier field A touches the maximum point of 
the holon dispersion, i.e., A c = 2tw c max |7(k)| = 6tw c . 
These collective charge excitations become gapped when 
A > A c , and a Mott insulating state appears, character- 
ized by (z a ) = with v = v = m = 0. 

Taking A = A c with v — v = m = 0, we can determine 
the quantum critical point from the following mean-field 
equations 
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Inserting w c from Eq. (J^U) into Eq. (|3"Tj) . one obtains 
the critical value for the interaction strength 
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= 0.312. (32) 



It is interesting to notice that the resulting param- 
agnetic Mott insulator has all kinds of lattice symme- 
tries. In particular, spin dynamics is described by gapless 
spinons. An effective field theory for spinon dynamics 
was proposed to be an SU(2) gauge theory with Dirac 
fermions.— 

It is not at all straightforward to understand dynamics 
of such gapless spinons due to complexity of the SU(2) 
gauge theory. It has been shown that an interacting sta- 
ble fixed point arises in the large- Nf limit ^ where Nf 
is the number of fermion flavors. Such a conformal in- 
variant fixed point was also shown to appear in the U(l) 
gauge theory with gapless Dirac fermions^ An interest- 
ing property of the stable fixed point is that the sym- 
metry of the original microscopic model, here the Hub- 
bard model, is enhanced, associated with special trans- 
formation properties of Dirac spinors i 14 ' 15 As a result, 
spin-spin correlations at an antiferromagnetic wave vec- 
tor have exactly the same power-law dependence as va- 
lence bond- valance bond correlations, which means that 
the scaling dimension of the staggered spin operator is 
the same as that of the valence bond operator^ This 
situation is completely unusual because scaling dimen- 
sions of these two operators cannot be the same in the 
level of the microscopic model. 

It is clear that one direct way to verify the algebraic 
spin liquid state is to observe the symmetry enhancement 
at low energies. If the staggered-spin correlation function 
turns out to display the same power-law behavior as the 
valence-bond correlation function, this will be an undis- 
putable evidence for the algebraic spin liquid phase be- 
tween the semi-metal phase and gapped spin liquid state. 
In the recent quantum Monte Carlo simulation data there 
seems to be uncertainty between the semi-metal phase 



and the gapped spin liquid state because such a simu- 
lation should be performed at finite temperatures. But, 
calculations for correlation functions need not be done 
at zero temperature. It is sufficient to show equivalent 
correlation behaviors in the quantum critical region at 
finite temperatures. 



B. From algebraic spin liquid to Z2 spin liquid 

Increasing more from the semi-metal to algebraic 
spin liquid critical point KcC7c , we find another param- 
agnetic Mott insulating phase, characterized by v 7^ 
and v 7^ with the d ± id' pairing symmetry. Re- 
call Eq. (27) and Eq. (28) for pairing symmetries 
that we checked explicitly. The algebraic spin liquid 
((z a ) = v = v = m = 0)to gapped spin liquid (v 7^ 
and v 7^ with (z a ) = m = 0) critical point is found 
with an ansatz of v = v — but v/v = 7^ 0. The 
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FIG. 2: (Color online) Spinon spectra ±D(k) given by Eq. 
(I25p with w = 1, v = 0.1 and m = 0. Left figure: If we 
consider only the d + id' pairing symmetry, spinon excitations 
are gapped only at the K point, but they remain gapless at 
the other K' point. Right figure: If we consider d + id' in one 
valley (K) and d — id' in the other valley (K'), spinon exci- 
tations become fully gapped. In this figure we assign d + id' 
to the upper band and d — id' to the lower band, respectively, 
in order to realize the d ± id' pairing symmetry. 
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mean-field equations to determine this critical point are 
given in appendix B. The system of equations is solved 
numerically, explicitly shown in appendix B. 

We find that the free energy reaches the lowest value 
for the d ± id' pairing symmetry. Actually, we checked 
self-consistency for various values of the angle parameter 
6 in Eqs. (27) and (28), and found 6 = ±tt/3 in Eq. (27), 
corresponding to d ± id'. In Fig. [2] we plot the spinon 
spectrum D(k) given by Eq. (|25|) for the d± id' pairing 
symmetry It shows that if only the d + id' pairing or- 
der parameter is taken into account, the energy spectrum 
opens a gap at one Brillouin zone edge (for instance, the 
K point), but it still keeps the Dirac cone at the other 
inequivalent Brillouin zone edge (the K' point). On the 
other hand, if we consider only the d — id' pairing sym- 
metry, we see that the K point remains gapless while 
only the K' point becomes gapped. This demonstration 
motivates us to assign the d + id' pairing symmetry to 
one valley (K) and the d — id' to the other (K 1 ), mak- 
ing the spinon spectrum fully gapped. Of course, this 
fully gapped state is energetically more favorable than 
the gapless state. In addition, this proposal resolves the 
problem of time reversal symmetry breaking at the same 
time. The edge state from the d+id' pairing in one valley 
is cancelled by that from the d — id' pairing in the other 
valley, preserving time reversal symmetry. One may re- 
gard this cancellation of such edge states as anomaly can- 
cellation due to fermion doubling in condensed matter 
physics. 

The critical value turns out to be K c U v /t = 0.315. Note 
that U v > U c . This intermediate phase between the semi- 
metal and gapped spin liquid is the algebraic spin liquid 
with an enhanced symmetry, as discussed in the previous 
subsection. We would like to emphasize that this region 
of U c < U < U v is not wide at zero temperature. But, 
the quantum critical region at finite temperatures will 
not be so narrow, and it will not be so difficult to ver- 
ify the algebraic spin liquid, considering staggered-spin 
correlations and valence-bond correlations. 

This pairing state can be verified by the quantized 
thermal valley Hall effecti 19 ' 20 The spinon number is not 
conserved due to particle-particle pairing, thus the charge 
Hall conductivity is not useful. On the other hand, both 
spin and energy (thermal) Hall coefficients arc important 
probes. But, the spin Hall conductivity vanishes due to 
the different assignment between two valleys. The ther- 
mal valley Hall effect should be observed in this state, 
regarded as the fingerprint of our Z2 spin liquid phase. 
This may be verified 2 ^ by either quantum Monte Carlo 
simulation^ or exact diagonalization 29 . 

It should be noted that our time reversal symmetry 
preserving Z2 spin liquid state is beyond the classification 
scheme based on the projective symmetry group because 
their possible Z2 spin liquids in the projective symmetry 
group are constrained with complete time reversal sym- 
metric pairingi^i In other words, d ± id' singlet pairing 
orders are excluded from the first although these pairing 
orders are not only found but also argued to be stable in 



recent studies! 1 ' 18 



C. From Z2 spin liquid to antiferromagnetic Mott 
insulator 

Our last subject is to investigate the quantum phase 
transition from the Z2 spin liquid to the antiferromag- 
netic Mott insulator. Here, we should take into account 
two order parameters such as the d ± id' pairing and 
antiferromagnetic ones. Generically, we expect four pos- 
sibilities. The first candidate is coexistence between such 
two orders, where the two critical lines cross each other. 
As a result, we have two critical points inside each phase. 
The second possibility is the multi-critical point, where 
the two critical points meet at one point. The third sit- 
uation will be the first order transition between them. 
The last corresponds to an intermediate state without 
any ordering, where the two critical points do not meet. 
First of all, we exclude the last possibility because this 
phase is nothing but the algebraic spin liquid and there 
is no reason for this reentrant behavior. 

We start to examine the possibility of coexistence. The 
antiferromagnetic critical point inside the Z2 spin liquid 
phase can be determined by m — while v and v are 
finite, thus determined self-consistently. The mean field 
equations for this quantum critical point are given in ap- 
pendix C-l. The strategy of solving the system of equa- 
tions is how to reduce the number of self-consistent equa- 
tions. Detailed calculations are provided in appendix 
C-l. As a result, we obtain two self-consistent equa- 
tions for two unknown variables. These equations can be 
solved numerically. For the first (n c = 1, k s = 1) and 
third (k c = 3/2, k s = 1/2) decomposition schemes in 
appendix A, we could show that there are no mean field 
solutions at the transition point. On the other hand, 
we find U m /t = 0.360 in the case of the d + id' pairing 
symmetry for the second decomposition scheme (k c — 1, 
k s = 1/2). 

The other quantum critical point is the Z2 spin liquid 
critical point inside the antiferromagnetic phase. It can 
be found when v ~ v — but m is finite, determined 
self-consistently. The mean field equations for this quan- 
tum critical point are given in appendix C-2. Solving the 
mean field equations self consistently, we could not find 
any solution. On the other hand, if the direct phase tran- 
sition from the antiferromagnetic Mott insulator to the 
semi-metal is concerned, we find the critical point occurs 
at U m /t — 0.330 for the second decomposition (k c = 1, 

Ks = 1/2). 

Our analysis for the quantum phase transition from 
the Z2 spin liquid to the antiferromagnetic Mott insu- 
lator shows that the nature of this transition depends 
on our phenomenological parameters of n c and k s . We 
could find the antiferromagnetic quantum critical point 
inside the Z2 spin liquid state for particular values of 
n c and k s while we could not obtain the Z2 spin liquid 
quantum critical point inside the antiferromagnetic Mott 
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insulating phase. We could not find the multi-critical 
point solution, either. As mentioned before, it is difficult 
to expect the algebraic spin liquid solution between the 
Z2 spin liquid and antiferromagnetic phases. Actually, 
we could find only one solution for the Z2 spin liquid 
to algebraic spin liquid transition, given by the previous 
subsection. The remaining possibility is the first order 
transition between the Z2 spin liquid and antiferromag- 
netic Mott insulator. We believe that the first order tran- 
sition is the generic case for the phase transition between 
these two phases. The formal procedure will be to inte- 
grate over spinons and holons and to obtain an effective 
Landau-Ginzburg- Wilson free energy functional for both 
d ± id' spin singlet pairing and antiferromagnetic order 
parameters. Based on the effective field theory, we can 
perform the renormalization group analysis and find the 
nature of the phase transition. This study is beyond the 
scope of the present study. 

One may ask the possibility of the Landau-Ginzburg- 
Wilson forbidden continuous transition between the Z2 
spin liquid with the d±id' pairing symmetry and the an- 
tiferromagnetic Mott insulator. Classification of Landau- 
Ginzburg- Wilson forbidden continuous transitions in two 
spatial dimensions has been performed in Ref. HO- Inves- 
tigating the classification table carefully, we can find that 
this transition does not belong to any cases. The main 
reason is that the singlet pairing order parameter cannot 
be symmetrically equivalent to the antiferromagnetic or- 
der parameter. The classification scheme reveals that 
the Neel order parameter can form a hyper-vector with a 
triplet pairing order parameter. In this respect we are al- 
lowed to exclude the possibility of the Landau-Ginzburg- 
Wilson forbidden continuous transition between the Z2 
spin liquid and the antiferromagnetic Mott insulator. 



IV. DISCUSSION AND SUMMARY 

In this paper we investigated the phase structure of the 
Hubbard model on the honeycomb lattice. Physics of one 
dimensional interacting electrons is our reference. As well 
known, even if we start from weak interactions, they be- 
come enhanced at low energies, destabilizing the Fermi 
liquid state. In one dimension such quantum corrections 
can be summed exactly, resorting to the Ward identity. 31 
The resulting electron Green's function shows two kinds 

SM ASL GSL AFM 

I • • * ► 

U c U v U m U/t 

FIG. 3: Schematic phase diagram. Abbreviations: SM is the 
semi-metal phase, ASL is the algebraic spin liquid, GSL is the 
gapped spin liquid, and AFM is antiferromagnetism. The SM- 
ASL and ASL-GSL quantum phase transitions belong to the 
second order while the GSL- AFM quantum phase transition 
is the first order. 



of branch cuts, corresponding to collective charge and 
spin excitations. In this diagrammatic approach it is 
difficult to see the nature of such fractionalized excita- 
tions. But, the bosonization approach is helpful at low 
energies, revealing that spinons and holons are identi- 
fied with topological solitons such as domain walls4 One 
can interpret this phenomenon in another respect that 
topological solitons acquire fermion quantum numbers 
via fermion zero modes, regarded as realization of quan- 
tum anomaly^ We believe that the spin-charge separa- 
tion in one dimensional interacting electrons results from 
not only just interaction effects but also hidden topolog- 
ical properties of Dirac fermions. Then, the next natural 
question is whether we can find this physics in higher 
dimensions. 

The graphene structure is an ideal system for realiza- 
tion of Dirac fermions. The first observation in this Dirac 
fermion system is that the vanishing density of states 
needs a finite value of the interaction strength U for an 
antiferromagnetic order to be achieved. Then, the ques- 
tion is whether we can find intermediate phases between 
the semi-metal and antiferromagnetic Mott insulator, al- 
lowing fractionalized excitations as one dimensional in- 
teracting electrons. Indeed, we could find two kinds of 
paramagnetic Mott insulating phases, which show frac- 
tionalized excitations. 

The algebraic spin liquid appears from the semi-metal 
state via the Higgs transition, gapping of charge fluc- 
tuations. Although it is not clear how the topological 
nature of Dirac fermions is introduced to result in such a 
spin liquid state, spinon excitations in the algebraic spin 
liquid can be identified with topological excitations cor- 
responding to meron (half skyrmion) excitations^ The 
underlying mechanism is that the symmetry of the origi- 
nal microscopic model is enhanced at low energies, allow- 
ing a topological term to assign a fermion quantum num- 
ber to such a topological excitation. The algebraic spin 
liquid turns out to have an 0(5) symmetry in the physi- 
cal case, where antiferromagentic correlations exhibit the 
same power-law dependence for distance as valence-bond 
correlations*^— It was pointed out that the correspond- 
ing effective field theory would be given by an 0(5) Wess- 
Zumino-Witten theory* 16 identifying spinons with such 
topological excitations. Comparing the algebraic spin liq- 
uid with the Tomonaga-Luttinger liquid, there is one to 
one correspondence between them except that charge ex- 
citations are critical in the Tomonaga-Luttinger liquid. 
Actually, spin dynamics in one dimension is governed by 
the 0(4) Wess-Zumino-Witten theory^ describing criti- 
cal dynamics of spinons. 

Because the stability of the algebraic spin liquid is not 
guaranteed beyond the large- Nf limit, we proposed how 
the quantum Monte Carlo simulation can prove the exis- 
tence of such a phase. As discussed before, the symmetry 
enhancement can be verified, calculating both antiferro- 
magnetic and valence-bond correlations at finite temper- 
atures. If such correlations turn out to have the same 
scaling behavior, we have the algebraic spin liquid phase 
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just beside the semi-metal state. 

When interactions are increased more, pairing corre- 
lations between nearest neighbor sites become enhanced 
in the singlet channel, destabilizing the algebraic spin 
liquid. As a result, spinon excitations are gapped due 
to their pairing orders. An interesting point is that the 
nature of this gapped spin liquid state is given by the 
d + id' singlet pairing order, which breaks time reversal 
symmetry. We would like to emphasize that time rever- 
sal symmetric combinations based on the d-wave pair- 
ing symmetry turn out to give higher energies than the 
d + id' pairing order. We suspect that this time reversal 
symmetry breaking may be related with the Berry phase 
effect of the momentum spaced One way to verify this 
statement is to check how the d + id' pairing symme- 
try is changed, increasing the chemical potential from 
the Dirac point, where the Berry phase effect becomes 
weaken. Unfortunately, the quantum Monte Carlo sim- 
ulation claimed that there is no time reversal symmetry 
breaking in the gapped spin liquid state. This incon- 
sistency was resolved, taking d — id' pairing to another 
valley. As a result, the edge state from the d + id' pairing 
is cancelled by that from the d — id' one. In addition to 
this time reversal symmetry, our proposal for the pairing 
order parameter turns out to be essential in order to have 
a fully gapped spectrum of spinon excitations, thus ener- 
getically more favorable than the case of only the d + id' 
pairing, where spin excitations remain gapless. We sug- 
gested an experimental signature, that is, the quantized 
thermal valley Hall effect as the fingerprint of this gapped 
spin liquid. 

Finally, we investigated the quantum phase transition 
from the Z2 spin liquid to the antiferromagnetic Mott 
insulator. We concluded that the first order transition 
will take place generically. We argue that this first order 
transition is involved with two symmetrically unrelated 
order parameters, displaying different discrete symmetry 
properties, here time reversal symmetry. We claim that 
the Landau- Ginzburg- Wilson forbidden continuous tran- 
sition will not appear, based on the existing classification 
scheme in the two dimensional Dirac theory on the hon- 
eycomb latticed 

We would like to point out that the SU(2) slave-rotor 
theory seems to overestimate quantum fluctuations. If 
one sets k c as the order of 1, the critical strength of the 
Mott transition is the order of 10" 1 for the critical value, 
compared with that from the quantum Monte Carlo sim- 
ulation. This overestimation originates from strong band 
renormalization for spinons and holons, given by effective 
hopping integrals, and Yy. Qualitatively the same 
situation also happens in the U(l) slave-rotor theory^ 
while the SU(2) slave-rotor theory seems to overestimate 
quantum fluctuations more. We believe that this aspect 
should be investigated more sincerely. 

Recently, the role of the spin-orbit interaction in the 
Hubbard model on the honeycomb lattice has been stud- 
ied both extensively and intensively, where one purpose 
is to reveal the interplay between the topological band 



structure given by the spin-orbit coupling and strong 
correlation effect. Novel exotic phases have been sug- 
gested in this Kane-Mele-Hubbard model, some of which 
are quantum spin Hall effect in a transition metal ox- 
ide such as Na2lr03^, a spin liquid state with a topo- 
logical band structure 3 -^, and the chiral spin liquid state 
with the anyon nature of excitations 36 . These interest- 
ing proposals will be verified based on " exact" numerical 
calculations 3 -!. 
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Appendix A: Decoupling scheme 

We discuss several decoupling schemes. The first ex- 
ample is 

i i 

+ fEE^-D' + fEE*-- 1 ) 

i <7 i a 

- f E( fjc t c -) 2 + |E n - ( A1 ) 

ia ia 

Formally, this magnetic decoupling does not correspond 
to the conventional Hartree-Fock analysis for antifcr- 
romagnetism because the interaction strength is twice 
larger the standard mean field value. 
The second possible decoupling is 

i i 

+ I EE d 2 + 1 EE ^-d 
+ \ EKE 4^)" - E «W]- (A2) 

i a a 

This decoupling recovers the standard mean field theory 
for antiferromagnetism, but the coefficient of the term 
XXE n io ~ I) 2 i s ViP 1 n °t equal to the coefficient of the 

i a 

term XXV^Vi) 2 - 

i 

The third possible decoupling is 

2H » = jE^^ + iE^l^) 2 

i i 

+ 7 EKE «W - E ™*W) 2 ]- (as) 

i a a 
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This equation determines A„. Once we find A„, we can obtain the critical value from Eq. (|B5|) together with Eq. 
(|Bip , given by 



3^El7(k) 
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Appendix C: To analyze the quantum phase transition from the Z2 spin liquid to the antiferromagnetic Mott 

insulator 

1. To find the antiferromagnetic quantum critical point inside the Z2 spin liquid state 

The mean field equations for this quantum critical point are given by 

.« Y^i-, i2_ 1 \ - w m | 7 (k)| 2 



,„ \^U|2_ / «c?7m W m \ ^ | 7 (k)| 2 ( 1 1 \ 

™V V 3 2iV/2^ J5(k) ^A m -2^(k) VA m + 2ti?(k)J' ^ j 



W|2_ / ^^m «m \ - [CWr f 1 1 \ 

t iV/2^ ^<| 7 (k)P+^|C(k)| 2 ' 



h ' ,U "' ' F( ; ' + , ± )• ( C6 ) 



3 JV/2 ^ - 2*£(k) V A ™ + 2 ^(k) 

Introducing x m = v m /w m and x m = v m /w m , we rewrite the above equations as 



\- ,2 1 \ - |7(k)|* 



4N/2 ^ ^| 7 (k)| 2 + ^|C(k)P' 



V^|A|2 = 1 \p ^m|C(k)| ,„ , 

s 4N / 2 u Vl7(k)| 2 + ^|C(k)P' 



K c C/ m 1 ^ | 7 (k)f 



I |2 / ^c u m r \ - 

« z. - v 6tfiJm 27V/2 V | 7(k) |2 + ^| C(k) |, 

1 



^-^(kTF+^iak)! 2 v /A m + v ^(kTF+5? 2 JC(k)P 



(C9) 
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f^c Urn 



E 



lC(k)| ; 



6tw m 2N/2 ^ ^| 7 (k)P +5^1^)1 = 



\ m - Vl7(k)| 2 +^|C(k)| 2 JX m + y/\>y(k)\2+a? m \ak)\ 



(CIO) 



^ K-sU m 1 \ 1 



(Cll) 



-. / k c U m 1 \ 



+ 



'X m - Vl7(k)| 2 +^,|C(k)| 2 V A ™ + Vl7(k)| 2 + ^|C(k)| 2 
where A r „ = 2tw m X m is redefined. From Eqs. (|C7[) and <|08J) we get 



£IC*I 5 

<5 



EN I 2 

5 



l y IC(k)| 2 

w / 2 k \/i^( k )i 2+2: ™ic( k )i 2 

1 l7(k)| 2 

w / 2 k /FKkTFSIJOkTF 



Similarly, Equations (fC9j) and (|C10|) give 



EIOI 2 



KM! 



W/2 k x/lT^P+^ICCk)! 2 \ v /A m - x /| 7 (k)P+£2 l |C(k)P v /A m + % /FKk)|2+£2j c(U )| 

lv72 E ' 



l7(k)| 2 



N/2 k x/lTC^P+^ICCk)! 2 \ v / A m -^|7(k)| 2 +5J l |C(k)| 2 v /A m +^/PKk)|2+£2 Ti | C ( k )| 

From Eqs. (jC9|) and ()C12j) we obtain 

l7(k)| 2 



1V72 E 



En 2 = ^ 



VlT( k )l 2+:r ™l^( k )l 2 V V /A ™-\ / lT( k )l 2 + 5; ™IC( k )l 2 A /A m + v ^(k)|2+2^|C(k)P 
JV72 E 



k V v /A ™-/FK k )l 2 +£^IC(k)| 2 v /A m + ^7k)?+5 2 „|C(k)|' 

Taking both sides of Eq. (|C9[) to the square power with w m from Eq. (|C7|) . we obtain 



■En 2 



K c U>n 



l7(k) 



iV/2 E ' 



l7(k) 



N ' 2 k Vl^WI'+^ICCk)! 2 

Inserting U m /w m from Eq. (|C11|) into this equation, we obtain 



tu. 



En 2 



l7(k)| 2 



N/2 k VlT( k )l 2+;r ™l?( k )l 2 I i/A^-^/RklP+S^ICtk)! 2 v /A m + ^7k)| 2 +£ 2 JC(k)| 2 
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Equations (|C15|) . (|C16j> . and (|C17j) give 
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(C13) 



(C14) 



(C15) 



(C16) 



(C17) 



(C18) 
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Inserting 



EN 2 ^E imi 



~ m i |7(k)|2 

^ vi7(k)i 2 +^ic(k)p 

from Eq. (|C13[) into Eqs. ()C14|) and (|C18|) . we obtain two self-consistent equations for two unknown variables, x m 
and X m . These equations can be solved numerically. We fix x m first, and solve these two equations for X m . Then, we 
obtain two functions, X m of x m . When two lines of these functions intersect, we obtain the solution of such equations. 
Once x m and X m are determined, the critical value of U m is also found from Eq. (|C11[) . 



2. To find the Z2 spin liquid quantum critical point inside the antiferromagnetic Mott insulator 

The mean held equations at this critical point are given by 



s AN ' 2 k ^/l7(k)| 2 + (^) 2 



*. E ICP = E W (~^— - ^== ) . (C2S) 



7(k) v v^w 



1 E / 1 ( C24 ) 



^?(7i=^ + 7rT^)' (C25) 

where A a — 2tw a X a is redefined. We solve these equations basically in the same way as the previous case. First, 
we reduce the system of six equations into two equations of two unknown variables, and solve the two equations 
numerically. 

From Eqs. (jC22|) and (fC25|) we get 



1 ^? l7(k)l yv^Mi: vzT^r }/ 

Wa = w 7 — (C26) 



5 N / 2 ? {vAw + vri^o) 

Equations (|C23)) and (|C25)) together with w Q from Eq. ([C20)) give 
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From Eqs. (|C21j) and (|C26I) we obtain 
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Equations (|C27|) and (|C28|) leads to 
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From Eqs. (|C24|) and (|C25I) we obtain 

6k s w a 



IV— 1 

^75 ^ yp^kjFRsp 



k V VAa-7(k) \Ao+ 7 (k) 



Inserting w a and u> a from Eqs. (IC20I) and (|C26|) into this equation, one obtains 

l7(k)| 2 1 v- I 



1 



1 y 

^T 5 k ^(k)FRS? 



1 V . 
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A/^a-7(k) V / Aa+ 7 (k) 7 ^ V V^" 7 (k) ^ V / Aa+ 7 (k) 



1 + 1 



(C29) 



(C30) 



(C31) 



Equations (|C29[) and (|C31| are the last two equations 
determining A Q and m a /tw a . We fix m a /tw a first, and 
solve the two equations for A a numerically. Then, we 



obtain two functions, A a of m a /tw a . When two lines of 
A a and m a /tw a intersect, we find the solution of these 
equations. 
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